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DISCONTINUOUS DUAL-PRIMAL MIXED FINITE ELEMENTS FOR ELLIPTIC 

PROBLEMS 


CARLO L. BOTTASSO*, STEFANO MICHELETTlt, AND RICCARDO SACCO* 

Abstract. We propose a novel discontinuous mixed finite element formulation for the solution of 
second-order elliptic problems. Fully discontinuous piecewise polynomial finite element spaces are used 
for the trial and test functions. The discontinuous nature of the test functions at the element interfaces 
allows to introduce new boundary unknowns that, on the one hand enforce the weak continuity of the 
trial functions, and on the other avoid the need to define a priori algorithmic fluxes as in standard 
discontinuous Galerkin methods. Static condensation is performed at the element level, leading to a 
solution procedure based on the sole interface unknowns. 

The resulting family of discontinuous dual-primal mixed finite element methods is presented in the 
one and two-dimensional cases. In the one-dimensional case, we show the equivalence of the method 
with implicit Runge-Kutta schemes of the collocation type exhibiting optimal behavior. Numerical 
experiments in one and two dimensions demonstrate the order accuracy of the new method, confirming 
the results of the analysis. 

Subject classification. Applied and Numerical Mathematics 

Key words, finite element method, mixed methods, discontinuous Galerkin, Petrov-Galerkin, 
elliptic problem 

1. Introduction and Motivation. We consider the following classical model problem: 

f — div(i/Vit) = / in 0 C R 2 , 
u = g on <9f2, 

where Dis a bounded polygon with Lipschitz boundary , = Oil, and u, /, g are given functions. 
Dirichlet boundary conditions are considered only for the sake of simplicity, but the extension to 
the case of more general boundary data is straightforward. Definitions of functional and geometrical 
entities are given in section 2. 

In view of the approximation of (1.1), we let {Th}h>o be a. family of triangulations of 0 and, for 
any 7ft, we denote by T any element thereof. Introducing the auxiliary unknown a = vXiu, problem 
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(1.1) can be reformulated in weak form over each element T of Th as 

/ v~ 1 a-Tdx+ / udivrdx — / ur-n T ds = 0, 

< Jj h' Jot ( 1 . 2 ) 

a • V vdx — v a • n T ds = fv dx, 

- Jt JdT JT 

where r and v are smooth vector and scalar functions, respectively, on T, while n T is the unit outward 
normal vector to dT. 

A unified framework for the analysis of the Discontinuous Galerkin (DG) method applied to elliptic 
problems in the form (1.2) was provided in ref. [1]. We shall follow the same path in the following, 
since the proposed method fits well in that framework. The DG formulation at the discrete level reads 
(see also [9]): find Uh G Uh and a h £ T h such that VT £ Th 


/ v 1 <Lh. ' Tft d x + / UhdivT h dx— Yj / h u T h -n T ds = 0 Vr ft £ S(T), 
J T J T _ AT " 6 


e€dT* 


/ Vb-Vvhdx- ^2 / v hk a ■ n T ds = / fv h 
* ' T „ Qrn *'6 


dx 


V« fc G t/(T), 


(1.3) 


eeeT-'s- 


where S(T) and G (T) are two sets of smooth vector and scalar functions, typically polynomials, defined 
on each element T, while and Uh are the corresponding global finite element spaces. We are now 
dealing with a Galerkin method where the values of the unknown fields on each edge e of dT , noted 
Uh \e and Q_ h \e, have been replaced by the numerical fluxes h u and h a . In practical implementations, 
h u and h a are chosen as functions of the internal unknowns Uh and a h on two neighboring triangles. 
In particular, it can be shown that most methods reported in the literature can be classified based on 
the particular choice of the numerical fluxes h u and h a , choice often inspired by ideas borrowed from 
finite volume methods. The interested reader can consult ref. [1] for exact details on the expressions 
of the fluxes in the different cases. 

It is clear that any particular choice on the nature of the numerical fluxes will affect the resulting 
scheme, eventually determining its numerical characteristics, stability, accuracy as well as the pattern 
of the stiffness matrix. This approach to the choice of the numerical fluxes is somewhat unsatisfactory, 
in the sense that it seems to be a rather artificial and dogmatic process, in some sense a “violence” 
to the underlying weak form (1.2). In this work, we shall therefore depart from the classical approach 
of defining a priori the expressions of the numerical fluxes. Hence, the interface fields h u and h a ■ n T 
are assumed to represent new additional problem unknowns , and we let the method implicitly define 
their values. In doing so, we obtain a genuine Petrov-Galerkin version of the DG method. 

This formulation generalizes the method described in refs. [2, 3, 4, 5] to multiple spatial dimen- 
sions. In the one-dimensional case, the method can be shown to be equivalent to collocation-type 
implicit Runge-Kutta (RK) schemes. In particular, when used in conjunction with Gauss-Legendre 
quadrature, the finite element method here proposed corresponds to the Kuntzmann-Butcher (Gauss) 
RK scheme, which enjoys well known properties of optimality. This provides a strong motivation for 
seeking an extension of those ideas to the multi-dimensional case, extension that is attempted in this 
work for the elliptic boundary-value problem. 

The paper is laid out as follows. Section 2 presents the proposed method for the classical problem 
of the Laplacian. Section 3 deals with the one-dimensional formulation of the present method and 
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with its connections with the time finite element and RK methods of refs. [2, 3, 4, 5]. In doing so, 
we shall temporarily abandon the simple linear model problem, and present the method for a general 
non-linear ODE problem. Section 4 is devoted to the description of the method in the two-dimensional 
case. Numerical results are discussed in section 5, where the proposed method is validated on one and 
two-dimensional test cases. Finally, the conclusions are discussed in section 6. 

2. The Discontinuous Dual-Primal Mixed Method for a Model Problem. We are now 

ready to define the proposed Discontinuous Dual-Primal Mixed (DDPM) discretization of the one- 
element weak formulation of problem ( 1 . 1 ): find Uh £ 1 /ft, & h € E ft , Aft £ Aft >9 and /Xft £ Aft such that 
VT £ 7ft we have 

o_ h ■ T h d . r + / Uftdivr^da;— / XhT h -n T ds = 0 Vr,. £ W(T). 

h j dT (2.i) 

/ CTft • Vt)ft dx - / S T VhHh ds = / f v h dx Mv h ^V(T). 

JT JdT Jt 

The meaning of the symbols is as follows. For any 7ft € {7ft}ft>o, let Eh be the associated set of 
edges. For any edge e of Eh, we denote by n e the unit normal vector to e, directed according to a 
uniquely fixed arbitrary convention. Namely, n e coincides with the unit outward normal vector on the 
boundary , if e fl , 7 ^ 0 , otherwise n e coincides with one of the two outward normals to the triangles 
sharing the edge e if efl , = 0. Then, for any T £ 7ft we denote by n T the unit outward normal vector 
to dT and define the sign function St = n e • n T = ± 1 . Through the use of St, we can associate the 
interface degrees of freedom with mesh edges, so that the sign function will determine whether a given 
flux Uh is entering of exiting from a given triangle. For simplicity, we shall omit in the following the 
dependence of n on the edge e which will be understood. 

Given the fact that the edge fluxes are now unknowns to the problem, we can not choose test and 
trial functions within the same spaces, as for the DG method (1.3). Therefore, we are now dealing 
with a Petrov-Galerkin formulation. The global finite element spaces for the internal variables q_ h and 
Uh are defined as 


Sft = {v h £ L 2 (f2) I v h \ T € S(T) VT £ Tft} , 

U h = {v h € L 2 (n) | v h \ T € U(T) VT £ r ft } , 

where Z/ 2 (fl) = (L^ifl)) 2 . Considering now the interface fields A ft and Uh, we note that they represent 
the traces of u and a - n T on dT, respectively. Using the language of computational mechanics, the 
new unknowns act as “mixed connectors” between neighboring triangles (see [16] and [18], Chapter 
13). We let A(e) indicate a set of scalar smooth functions defined on each edge e. Then, we define the 
(global) finite dimensional spaces 

A ft = {rth £ L 2 (£ h ) \j] h \g £ A(e) Ve £ £ h } , 

Aft, ff = {Vh £ Aft | % = if e n , ^0}, 

where P is the T 2 -projection operator on A(e). Finally, we introduce the local and global finite element 
spaces for the test functions r and v. Denoting by W_(T) and V (T) two sets of smooth vector and 
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scalar functions defined on each element T, we have 


W h = {r h 6 L 2 ( 12) | r h | T € H^T) VT 6 T ft } , 

^ = {w h e T 2 (12) I Wftlr e E(T) VT e T h } . 

The choice of the functional spaces will be detailed in the following sections for the one and the 
two-dimensional case. 

We note that the method is completely conservative in the sense of ref. [1], since the interface 
unknowns are associated with edges and are consequently identical for two triangles sharing the edge. 
This last property implies the following conservation statement 

Phds + I f dx = 0, 

where U is the union of some collection of elements, and Vh was taken to be identically unity in the 
second of (2.1). 

The DDPM approach can be viewed as a consistent hybridization of the dual-primal mixed (DPM) 
method introduced in ref. [12]. Indeed, the spirit of the standard hybridization procedure (see [7], 
Chapter V) is to relax the continuity of some of the unknowns (a • n in the dual hybrid philosophy or 
u in the primal hybrid philosophy) and then to enforce it back through the introduction of suitable 
inter-element Lagrange multipliers. Precisely, u\ L represents the Lagrange multiplier in the case of 
dual hybrid methods, while a • n\e is the Lagrange multiplier in the case of primal hybrid methods. 
In the present case, we relax the continuity of both u and a ■ n in the interior fields and let the 
interface variables connect neighboring elements. The interface variables A and p do not play the 
role of Lagrange multipliers in our case but, nevertheless, enjoy appealing properties. Physically, they 
have the clear meaning of u\e and a ■ n\ L , respectively, while numerically they exhibit a higher rate of 
convergence than the internal variables, as it will be shown in the following. This is also a feature of 
the Lagrange multipliers in the case of standard hybrid methods ([7], Chapter V). 

Alternatively, one could note that, by gluing the test functions together at neighboring triangles, 
each internal edge receives equal and opposite contributions from the two elements using it, and 
therefore the boundary unknowns are eliminated. This way one recovers from (2.1) the DPM method. 

3. The DDPM Method in the One-Dimensional Case. We abandon now for the moment 
the linear model problem (1.1), and we consider a generic system of first order ODE’s 



z' = f(z,x), (3.1) 

where z,f£ R U and D is the number of state variables. Suitable initial or two-point boundary 
conditions complement problem (3.1). In the case of the linear model problem, we clearly have 
z = ( u,a) T and / = (cr, — f) T . However, since the discussion can be easily completed for the more 
general case, we consider (3.1) instead of the one-dimensional version of (1.1) in the remainder of this 
section. 

We assume 12 = (0, A') and let 0 = xo < X\ < ... < x n -i < x n = X be a given partition Th of 
12 into n > 1 intervals T) = [xi,Xi + 1 ] of size hi, i = 0, . . . ,n — 1, where {xi}f =0 are the nodes. Note 
that in the one-dimensional case £/, = {xi}f =0 and each edge e degenerates into one single node. 
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The statement of the present finite element method is as follows: find (z h ,X h ) G (Z h x A ft ) such 
that for each T G 7h, « = 0, . . . , n — 1, the following holds 



• w h dx - (A i+1 


■w h {x i+1 ) - \i-Wbixi)) =0 


Vw h eW{Ti). (3.2) 


For any T &Th, we define the local finite element spaces for the internal fields z h and for the test 
functions as follows: 


Z(T) = (F k (T)) D , W(T) = (P k+1 (T )) D , 
with k > 0. The local finite element space for the interface unknowns \ h is 

A(3Tj) = {Aj, A i+1 }, t = 0, . . . , n — 1, 

where for any function A G A(9Tj) we let Aj = A(#i) and A*+i = A(®*+i)> * = 0, . . . ,n— 1. Furthermore, 
we set 


U(T) = U(T) x A (3T), V(T) = W(T). 

The finite element space is then 

OP k (T) = {U, A; w) | U,A;w) e U{T) x V(T) VT e T h }. (3.3) 

The global finite element spaces for the internal fields and the test functions are 

= {z h € (L*m D | z h | T € Z_(T) VT €%}, 

W h = {w h G ( L 2 (Q)) D | w h |t G W_(T) VT G Th) , 

while the global spaces for the interface unknowns are 

A fc = {{A,}2 =o |A,eR D ,t = 0,...,n} ) 

where functions are defined only at the nodes of Eh- Finally, gathering the previous definitions, we 
obtain the global finite element space as 

OP£ = (Z h X A h ) x W h . (3.4) 

It is clear that both initial and boundary-value problems can be solved within this framework. 
For an initial value problem, one can solve marching sequentially in an element-by-element fashion, as 
with any time stepping procedure. We have: dim(Z(T)) = D(k+ 1), dim(W(T)) = D(k + 2) VT G Th- 
Therefore, for any T G Th, (3.2) leads to a system of D(k + 2) equations in the D(k + 1) + 2D 
local unknowns z h [r, X h \x, leaving space for the D initial conditions that complement problem (3.1). 
A boundary-value problem leads to a global discrete problem defined over f l. In this case we have 
nD{k- 1-2) equations in the nD(k+l)+2nD unknowns. However, considering that A,+i appears in both 
neighboring elements T and T + i , the total unknown count is reduced to nD(k + 1) + 2nD — (n — 1)£>, 
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which leaves once again space for D two-point boundary values. In both cases, the internal field z h 
can be eliminated at the element level in favor of the interface unknowns A ft , and then recovered at 
convergence. 

We note that the scheme presents two jump discontinuities for each finite element, namely one at 
X{ and the other at Xj+i, since A* ^ z h (xi ) and A«+i ^ z h (xi+i). For this reason it was termed the 
Bi-Discontinuous (BD) scheme in ref. [3]. It is clear that such symmetric treatment of the element 
boundaries implies no privilege in the direction of the flow of information, and it is therefore best 
suited for two-point boundary-value problems. 

Quite differently, for initial value problems it is sometimes desirable to have stiffly accurate 
schemes, that are able to damp the higher frequencies components from the computed response [11]. 
This is achieved in practice by forcing a lack of symmetry in the scheme, that incorporates the knowl- 
edge on the direction of flow of information within the element. It is possible to derive such schemes 
in the present framework, allowing one single jump discontinuity at Xi while enforcing the condition 
Aj+i = z h (xi. (- 1 ). One obtains in this case the so called Singly-Discontinuous (SD) schemes [3]. 

3.1. Equivalence with Runge-Kutta Methods. It was shown in ref. [3] and then again with 
a slightly different proof in ref. [4], that the finite element method (3.2) can be written as a RK process 
for any order k. RK methods are probably best known for the solution of initial value problems, but 
clearly can also be used for the solution of two-point boundary-value problems by assembling the 
single steps to yield a global discrete problem defined over the whole computational domain, exactly 
as in the finite element method. Since the link between the two approaches seems to be useful in the 
characterization of the proposed method, we briefly review this result in the following. 

We begin by selecting a quadrature rule for the evaluation of the integrals in (3.2). The rule is 
defined by s abscissae c, and weights 6*. We consider the case s = k + 1, therefore the number of 
quadrature points is the same as the number of finite element nodes of the trial functions. 

Following the usual FEM practice, one expresses the discrete equations in terms of the nodal 
values. Here we shall take a different approach, and assume as unknowns the values of the finite 
element trial functions at the quadrature points. This is the key idea, for showing the link existing 
between finite element and RK methods. The change of unknowns can be done because we are using 
a quadrature rule that uses as many integration points as the number of finite element nodes, as 
previously said. We will come back later to some interesting consequences of this assumption. 

We first note that the test functions can be expressed as 


k + 2 

= E - < " 1 

i= 1 


QLi, 


(3.5) 


where are polynomials, r € I, and a k their associated amplitudes. Next, let us define the (s + 1) x s 
matrix P, 'c = [r* -1 | r= ^] , with i = 1, . . . , s + 1, j = 1, . . . , s, the fj’s being s real values. The derivative 
of with respect to r is then P £ = [(j — l)r* _2 | T= ^.] . Note that, for any choice of the fj’s, the first 
row of P , f is all made of one’s, and the first row of P ^ is all made of zero’s. For fj = cj, we have the 
matrices P c = [cf _1 ] , P' c = [(* — l)c*- — 2 ] . We also define the s x s diagonal matrix B = N , with 
* = 1, . . . , s, and the following .s-dimensional vectors: b = (&i, & 2 , . . . , b s ) T , c= (ci,c 2 , . . . , c s ) T , and 
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(3.6) 


0 = (0, 0, , 0) T , 1 = (1, 1, ... , 1) T , Ij = (1,0,... , 0) T . 

We note that 


P'B = 


r 

Q' 


Q being the s x s matrix Q = [r*| T=Cj bj\ = \cjbj] . The derivative of Q with respect to r is given by 
Q' = \ic l f~ l bj\. Furthermore, if the quadrature formula is of order s, we have that 


Q'l = I- 


(3.7) 


Equation (3.7) corresponds to the so called “B simplifying assumption'’ in the theory of RK methods. 
The discrete weak form (3.2) can now be written as 


{{P' C B) ®I D )z + hi((P c B) <g> I D ) l = (1 ® Id) A i+1 - (^ <B> In) A*, 


(3-8) 


where the following two sD-dimensional vectors were defined: z = {zf,-.. ,1®) T , / = (/(l 1 ,®* + 
Cihi), . . . , f(z s ,Xi + c s hi)) T . Id is the D x D identity matrix, while <g> denotes the tensor Kronecker 
product of matrices, so that (•) <g> I# ensures that all degrees of freedom of the vectorial problem are 
integrated according to the same rule. 

Solving for z and A i+1 and using (3.7), we have 


(l,A j+1 ) T = 


1 hi(lb r -Q r Q) 

1 hib T 


'MA i,f) T . 


(3.9) 


From this expression, we conclude that (3.2) can be written as a s-stage RK method whose tableau 
[10] is given by 


1 b T -Q'-'Q 


(3.10) 


Finally, we test the Gauss-Legendre, Lobatto and Radau quadrature rules and compute the cor- 
responding tableaux (3.10). The results are summarized in Table (3.1). 

Table 3.1 

Correspondence between quadrature rules for BD finite elements and RK methods. 


Quadrature Rule 

RK Method 

Gauss 

Kuntzmann-Butcher 

Lobatto 

Lobatto IIIB 

Radau-Left 

Radau IA 


Table (3.1) shows that the family of schemes deriving from the one-dimensional finite element 
method considered in this work corresponds to some well known RK algorithms. The two approaches 
differ on the choice of the unknowns: the finite element method solves for the nodal values , while 
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the RK method solves for the quadrature point values. However, these two sets of values are simply 
related by a linear transformation. Furthermore, note that all RK schemes appearing in the table 
are of the collocation type, where polynomials interpolate the two interface values \ i , A i+1 and the 
internal stages z [10] . This link between finite elements and RK processes allows a somewhat unique 
way of looking at the discretization procedure: from the finite element point of view we see jumps 
in the field variables, since we have discontinuous interpolations of the internal fields, glued together 
from one element to its neighbor through the presence of the interface unknowns. Quite differently, 
from the RK point of view we have continuous interpolations of the interface unknowns with internal 
variables at the collocation points, with apparently no jump discontinuities in the solution. It should 
also be remarked that for these RK methods, maximal order is in general attained only at the interface 
values, and not at the internal stages. We then expect to observe higher rates of convergence of the 
interface fields also in the multi-dimensional generalization of the method. 

The same table shows that the present finite element formulation, when used in conjunction 
with Gauss-Legendre quadrature, yields the Kuntzmann-Butcher (or Gauss) RK methods. These RK 
schemes are probably the ones that enjoy the greatest number of numerical properties, and are in 
this sense optimal. In fact they are of maximal order (2s), they are algebraically stable, and also 
symplectic when applied to Hamiltonian systems [17]. Simplecticity is the distinguishing property of 
Hamiltonian systems, and its numerical preservation has a strong influence on the behavior of the 
integration scheme. 

As a final point, we note that these results are solely based on having assumed a quadrature rule 
that uses as many points as the number of finite element nodes. In reality, working from the point 
of view of the finite element method, one could choose to use a greater number of quadrature points, 
for example in the hope of improving the accuracy in the evaluation of strongly non-linear functions 
in the term / f(z h ,x) ■ w h dx. It can be shown however, that increasing the number of quadrature 

JTi ~ 

points can in reality degrade the performance of the method. For example, using a larger number of 
points destroys the symplectic nature of the scheme [3]. 

4. The Two-Dimensional Case. For two spatial dimensions, the choice of the finite element 
spaces for the internal and interface unknowns is a straightforward extension of those introduced in 
the one-dimensional case. Precisely, for k > 0 we set 

S(T) = (P ft (T)) 2 , U(T) = F k (T) VT&T h , 

A(e) = Pfc(e) Vee£ ft . 

The choice of the finite element test spaces is less obvious and requires some care. Here we 
focus on the description of the method in the cases k = 0 and k = 1. We start by introducing the 
Raviart-Thomas finite element space of degree k [14] 

ET * (T) = (P* (T )) 2 © *P* (T) VT G Th- 
in the case of the OP° finite elements, we set W(T) = MT’o(T') and V(T) = Pi(T). The global finite 
element spaces W h and V k are then simply defined as the products of the corresponding local spaces. 
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Notice that W_ h and Vft are the discontinuous counterparts of the corresponding test spaces that have 
been used in the dual-primal formulation of ref. [12]. 

In the case of the BP^ finite elements, we set W_(T) = KTi(T) and V(T) = ^(T) © B 3 (T) for 
each T G 7ft. B 3 (T) is the following cubic bubble [15] 

B 3 = (<Po ~ <Pi)(<Pi ~ ‘Pl)( l P 2 ~ 7>o), 

where tpi = i = 0, 1,2, are the centroidal coordinates of a point iel ! with respect to the 

vertices of T. Enrichment of the scalar test space is a standard procedure in hybrid methods when 
the degree of the polynomials in V (T) is even (see refs. [15, 8]). Actually, one can check that without 
adding the cubic bubble the rectangular matrix arising from the integral / Sr i'h l->h ds becomes rank 

JOT 

deficient. 

As far as the implementation is concerned, we point out that for both the k = 0 and k = 1 cases 
it is possible to perform a static condensation of the internal variables in favor of the edge unknowns, 
obtaining a linear system in the sole interface variables Aft and //ft. 

The generalization of the DDPM method to higher orders will be more extensively addressed in a 
forthcoming paper. Let us just mention that a possible general recipe for constructing the DP^ finite 
element spaces for even k can be given as 

W(T) = BDFM k+1 (T), V(T) = P fc+1 (T) VT G T h , 

where BDFM k + 1 is the triangular analogue of the reduced element introduced in ref. [6]. 

£ 

5. Numerical Results. In this section we demonstrate the OP ft finite elements for k — 0 and 
k = 1 through the solution of test cases in both one and two spatial dimensions. 


5.1. One-Dimensional Examples. We consider the following two-point boundary-value model 
problem 


-(i '(x)u'(x))’ = f(x), 0 < x < 1, 

u(0) = u(l) =0. 


In all numerical experiments, successively finer uniform grids of size h = 2“ J , j G {1, ... ,9} are used. 
We set e u = u — Uh and e a = a — <j h ; moreover e\ and denote the functions defined only at the 
nodes of E k such that e\(xi) = u(xi ) — and e^Xi) = a (a^) —//,,* = 0, . . . ,n. 

For any function v G L 2 (fl) and sufficiently smooth on each T £ Th, we denote by I] //’ft’) the 
approximation of f T v dx using the trapezoidal rule. Moreover, for each T G 7ft we denote by {Xj T }j=o 
the k + 1 Gauss-Legendre points on T. In order to measure the convergence rate of the DP^ finite 
elements, we use the following three norms: 

( n — 1 
i = 0 



|M’||ft,GP,oo = 


max I«;(a:n7’))|, 
ie{ 0 , •••,»-!} u ’ 


max 


»e{o, •••,»-!} \je{o,i} 


max (\w(x j<T )\) , 


if k = 0, 
if k = 1, 


IMIft,oo= max 

,n} 
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where w is any bounded function defined at the Gauss points of Th and r] is any bounded function 
defined at the nodes of Th- Clearly, || • ||fc j0 ,n is a discrete L 2 (Cl) norm, while the other two norms are 
discrete maximum norms over the set of Gauss points and £h, respectively. The quadrature formula 
used in the discrete norm is accurate enough not to pollute the computed accuracy, and at the same 
time avoids to sample the solution at super convergent points. Moreover, observe that any integral 
appearing in the DDPM discretization of (5.1) has been computed using a one-point Gauss-Legendre 
quadrature rule (with node ,r[] on each T e Th) for k = 0 and a two-point Gauss-Legendre quadrature 
rule (with nodes {art t}}=o on each T e Th) for k = 1. 

The following symbols are used in the graphs: for any h , symbol refers to ||e M ||h,o,fi (or 
||e CT |U,o,n), while symbols and ‘o’ denote ||e u ||A,Gi’,oo (or IKHa.gp.oo) and ||e A || fti00 (or ||e M || ft>00 ), 
respectively. 


Test Case 1 We set v{x ) = 1 and f(x) = sin(7nc), x € [0,1]. The exact solution is the smooth 
function u(x) = (1/7T 2 ) sin(7rai), which implies a(x) = u'(x) = (1 /tt) cos(irx). We show in figure 5.1 
the error curves for e„ and eA (left) and e a and (right) in the case k = 0. The corresponding error 
curves in the case k = 1 are shown in figure 5.2. 




Fig. 5.1. One spatial dimensions, test case 1. Approximation errors e u and e\ (left) and e a and e ^ (right) for 
DP h finite elements. 

The plots show second-order convergence for the interface unknowns in the case k = 0, while the 
internal unknowns are only first-order accurate. Superconvergence of the internal variables can be 
observed at the Gauss point. For the DPjj finite elements, the error curves show that the interface 
unknowns are fourth-order convergent, the internal fields are second-order accurate while a third-order 
convergence rate is exhibited by the internal unknowns at the two Gauss points of each element. All 
these results are in agreement with the theory of RK methods. 
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Fig. 5.2. One spatial dimensions, test case 1. Approximation errors e u and e\ (left) and e a and e y (right) for 
OP /j finite elements. 


Test Case 2 This test problem demonstrates the performances of the DDPM method in the presence 
of strong solution gradients. The coefficient v and the source term / are chosen as in ref. [13] (see 
also ref. [8], section 4.4.1, for further comments), and read 

v(x) — — + a( x - x) 2 , 
a 

f(x) = 2 (1 + a(x - x) {tan _1 (a(x - x)) + tan~ x (ax)}) . 

The exact solution is 

u(x) = (1 — x) (tan -1 (a(x - x)) + tan -1 (ax )) , 

a(x) = ( — f- a(x — x) 2 ) { — tan -1 ( a(x — x)) — tan -1 (ax) + (1 — x) zrnr 1 , 

\a ) \ l + a z (x — x) z J 

where a(x) = v(x)u'(x). Furthermore, we set a = 100 and x = 0.36388. 

We show in figure 5.3 the error curves for e u and e\ (left) and e a and e M (right) in the case k = 0. 
The corresponding error curves in the case k = 1 are shown in figure 5.4. 

The numerical results clearly show the increased “roughness” of the problem when compared 
with the corresponding error curves for test case 1. Nevertheless, all variables attain the expected 
convergence rate as h — )■ 0. Figure 5.5 shows the exact solution u (solid line) plotted against the finite 
element solution computed with OP^ elements and h = 1/20. Notice the sharp internal layer, 
resolved within one element. 

5.2. Two-Dimensional Examples. For the two-dimensional case, we consider (1.1) in the 
square domain fl = (—1, l) 2 with g = 0. The exact solution is in this case u(x, y) = (x 2 — 1 )(y 2 — 1) and 
(?(x,y) = 2 (x(y 2 — l),y(x 2 — 1)) T ; the right-hand side / is computed accordingly. In all numerical 
experiments, we use a uniform grid of isosceles right triangles of size h x = h y = h , where h takes on 
the values 2 1_ - 7 /5, j € {0, 1,2,3} for the JW° h case and {2/5,2/10,2/20,2/30} for the OP), case. 


li 






Fig. 5.3. One spatial dimensions, test case 2. Approximation errors e u and e\ (left) and e a and e ^ (right) for 
\ finite elements. 



Fig. 5.4. One spatial dimensions, test case 2. Approximation errors e u and e\ (left) and e<j and e ^ (right) for 
OP h finite elements. 


As in the one-dimensional examples, we denote by e u = u — Uh and e a = a — a h the discretization 
errors associated with the internal fields, while e\ = V\ — A/, and e M = P/x — Hh are the discretization 
errors associated with the interface fields. V is the L 2 -projection operator on A(e), for each e e £h- 
Furthermore, for any T G Th we denote by |T| the measure of T, by a r , r = 0, 1, 2 the vertices of T, 
by x C T the centroid of T and by x^ T ,j = 0, 1, 2 the three Gauss points of T obtained from the area 
coordinates {2/3, 1/6, 1/6} by permutation. In order to measure the convergence rate of the OP£ 








Fig. 5.5. One spatial dimensions, test case 2. Exact solution (solid line) and numerical solution A/j (‘o’) computed 
using MP ^ finite elements (h = 1/20 ). 


finite elements, we introduce the following discrete norms: 


IMko.fi — 

IMk GP,oo 




KTeTh 


r = 0 


m&x\v(x CT )\, 
max max Jv(^ T )|, 


T&Thj€{ 0 , 1 , 2 } 


1/2 


\\v\\h-l/2,£ h 


Ek 

e€Sh 


if k = 0, 
if k = 1, 


where v is any sufficiently smooth function on T belonging to L 2 (Q), r] is any function belonging to 
I/ 2 (e) for each cefj and |M|o,e the Z/ 2 -norm on edge e. Notice that || • lko.fi is a discrete I/ 2 (f2) 
norm, || • \\h,GP,oo is a discrete L°°(f2) norm, while || • lk-i/ 2 .£,, is a discrete i? -1 / 2 (£ft) norm (see also 
ref. [7], Sect. V.3). For vector functions v = (i>i,f 2 ) T , the above norms are defined 

( A 1 / 2 

Nko.fi = ^Xu=i,2 Nlko.fi J > 

Nllft.GP.oo = max i6 { li2 } Nllft.GP.oo- 


In the following figures, symbol l *’ is used to denote ||e„|ko,fi (° r lle^Uft.o.o), while symbols 
and ‘o’ denote ||e„|k G p,oo ( or INIU.gp.oo) and ||e A |U,_i/ 2 ,£ h (or ||e M || ft _ 1/2)£h ), respectively. 

We show in figure 5.6 the error curves for e u and e A (left) and e a and e M (right) in their respective 
norms for the case k = 0. The corresponding error curves in the case k = 1 are shown in figure 5.7. 

For the u field, the plots show second-order convergence for the interface unknowns in the case 
k = 0, while the internal unknowns are only first-order accurate. Superconvergence of the internal 
variables can be observed at the Gauss point. This behavior is similar to the one observed in the 
one-dimensional case. For the a field, the plots show once again second-order convergence for the 
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Fig. 5.6. Two spatial dimensions. Approximation 
elements. 



Fig. 5.7. Two spatial dimensions. Approximation 
elements. 



errors e u and e\ (left) and e g and e ^ (right) for DP h finite 



errors e u and e\ (left) and e a and e ^ (right) for DP h finite 


interface unknowns, while the internal unknowns are only first-order accurate. However, for this field 
superconvergence at the Gauss point is not observed. 

In the case of DP^ finite elements, the error curves for the it field show third order accuracy 
for the interface unknowns, and second order accuracy for the internal unknowns. For the a field, 
the interface unknowns achieve an order equal to approximately 5/2, while the internal variables are 
second-order accurate. Once again, superconvergence is not observed at the Gauss points. 
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In all cases, the interface unknowns are consistently more accurate than the internal fields, as 
expected from the one-dimensional analysis, although the exact order achieved by the various fields 
does not seem to obey an obvious rule. As far as the order of convergence of Aft is concerned, the 
results agree with the classical convergence analysis of dual- hybrid methods. Actually in this case 
one obtains 0(h k+ 2 ) when the Raviart-Thomas finite elements of degree k are employed (see [7], sect. 
V.3, formula (3.19)). 

Let us conclude by comparing the performance of the DDPM approach to the performance of the 
standard dual method. In particular, we consider the DP° method with static condensation versus the 
dual method using lowest order Raviart-Thomas finite elements (D°) with no hybridization. Table 5.1 
shows the results obtained for the example here considered. The meaning of the various entries is the 
following: h is the value of the mesh size on the boundary of the domain, dim is the dimension of the 
matrix associated with the linear system, nnz is the number of non-zero elements of the same matrix 
and flops is the number of floating-point operations as computed by Matlab using sparse operations. 
The table on the left refers to the DP° case, while the table on the right collects the results for the 
0° h case. In the first case, flops includes assembly, solution of the linear system and recovery of the 
internal fields, while in the second case it considers the sole assembly and solution phases, since there 
are no recoveries to perform. 


Table 5.1 

Flop counts for the DP h and D ft methods. 


h 

dim 

nnz 

flops 

h 

dim 

nnz 

flops 

2/5 

150 

662 

41424 

2/5 

135 

617 

26888 

2/10 

600 

2713 

224775 

2/10 

520 

2272 

209964 

2/20 

2400 

11060 

1835876 

2/20 

2040 

9752 

2993179 

2/40 

9600 

44267 

16867719 

2/40 

8080 

36112 

22407438 


The tables show that the DP° method performs better than the D° method in terms of flop-count 
as the dimension of the linear system grows, despite the fact that the dimension of the matrix and 
the number of its non-zero entries are greater than in the D° case. 

6. Conclusions. We have presented a novel mixed dual-primal finite element formulation for 
elliptic boundary-value problems. This work is a first attempt at generalizing discontinuous finite 
element formulations for generic ODE problems that were developed in previous papers. The one- 
dimensional formulation can be shown to lead to optimal RK methods, and therefore provides a strong 
incentive towards the generalization to multiple space dimensions. 

The DDPM method uses interpolations of the unknown fields which are discontinuous between 
neighboring triangles, in the same spirit as the DG method. However, no particular expression of the 
numerical fluxes is selected a priori in this case. Additional unknowns are introduced at the element 
interfaces, that act as connectors gluing neighboring triangles together. The discontinuous internal 
unknowns can be eliminated at the element level, leaving a solution scheme in the sole interface 
variables. A higher rate of convergence is observed for the interface unknowns with respect to the 
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internal variables, as expected from the analysis in the one-dimensional case. Numerical examples 
have been presented to support the analysis and provide numerical evidence of the characteristics of 
the method here investigated. 

Future efforts will concentrate on the extension of these ideas to other model problems, and on 
comparisons of this class of methods with alternative, well established solution procedures. 
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